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We studied how the period length and the mass ratio affect the thermal conductivity of isotopic nanoscale 
three-dimensional (3D) phononic crystal of Si. Simulation results by equilibrium molecular dynamics show 
isotopic nanoscale 3D phononic crystals can significantly reduce the thermal conductivity of bulk Si at high 
temperature (1000 K), which leads to a larger ZT than unity. The thermal conductivity decreases as the 
period length and mass ratio increases. The phonon dispersion curves show an obvious decrease of group 
velocities in 3D phononic crystals. The phonon's localization and band gap is also clearly observed in spectra 
of normalized inverse participation ratio in nanoscale 3D phononic crystal. 



Thermoelectric materials are important for generating electricity from waste heat and being used as solid-state 
Peltier coolers. The performance of thermoelectric materials depends on the figure of merit (ZT) 1 , ZT = 
S 2 aT/K, where S, T, a, and k are the Seebeck coefficient, absolute temperature, electrical conductivity and 
total thermal conductivity, respectively. ZT can be increased by increasing S or a, or decreasing k. However, it is 
difficult to improve ZT in conventional materials. First, simple increase S for general materials will lead to a 
simultaneous decrease in a 1,2 . Also, an increase in a leads to a comparable increase in the electronic contribution 
to K 1,2 . 

An alternative way to increase ZT is to reduce the thermal conductivity without affecting electronic property 1,3 . 
Moreover, ultra-low thermal conductivity is also required to prevent the back- flow of heat from hot end to cool 
end. Therefore, reduction of thermal conductivity is crucial in thermoelectric application. 

Phononic crystals are constructed by a periodic array of scattering inclusions distributed in a host material. Due 
to its periodic change of the density and/or elastic constants, phononic crystals exhibit phononic band gaps 4 . This 
remarkable property is very different from those of the conventional materials and can be engineered to achieve 
new functionalities. A special one- dimensional phononic crystal is the superlattice, one dimensional periodic 
arrangement of two different materials. It is demonstrated that superlattice crystals are effective to achieve very 
low thermal conductivity 5 " 11 . Superlattices have been extensively studied to design thermoelectric materials with 
high ZT. Preliminary works show that there is a minimum value of thermal conductivity in the direction 
perpendicular to the planes of superlattice when the period length is reduced to nanoscale 12 " 14 . 

Recently, it is demonstrated experimentally that the Si nanomesh film, a two-dimensional (2D) phononic 
crystal, exhibited low thermal conductivity 15 by modification of phonon band structure. Single crystalline Si by 
phononic crystal patterning in 2D has a smaller value of thermal conductivity (~6 W/m-K) than bulk Si because 
of the low group velocities and the coherent phononic effects 16 . It is predicted that atomic-scale 3D phononic 
crystal of Ge quantum- dot in Si has very low thermal conductivity in all three spatial directions 17 . The thermal 
conductivity is reduced by several orders of magnitude compared with bulk Si. This reduction of thermal 
conductivity is due to the reduction in group velocities and multiple scattering of particle-like phonons. 

In this letter, we study the thermal conductivity of nanoscale 3D silicon phononic crystal. The 3D crystal 
consists of 28 Si atoms and "isotopes" M Si atoms which have the same properties as 28 Si except the mass, where M is 
the atomic mass of the isotope of Si. The mass ratio, R, is defined as R = M/28. The 3D isotopic phononic crystal 
could also be named as 3D superlattice because different material arranged periodically in three spatial directions. 

We find the 3D isotopic phononic crystal has the ability to flatten phonon dispersion curves compared with 
that of bulk Si and it could show band gaps when properly arranged. We studied how the period length and the 
mass ratio affect the thermal conductivity of the 3D phononic crystal. The phonon dispersion curves and inverse 
participation ratio are also computed to understand the mechanism of the reduction of thermal conductivity. 
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On the other hand, the scatterings of isotopic doping could sig- 
nificantly reduce the lattice thermal conductivity without much 
reduction of the electrical conductivity 114 . As a result, the thermal 
conductivity of 3D isotopic phononic crystal has a quite low value, 
which can lead to a larger ZT than unity. 

Results 

Fig. 1 shows the structures of the isotopic nanoscale 3D phononic 
crystal, where 28 Si and M Si atoms are assembled periodically in three 
spatial directions. Fig. l(a)-(d) shows the structures of 3D phononic 
crystals with different period lengths, corresponding to 1.09 nm, 
2.17 nm, 3.26 nm and 6.52 nm, respectively. The volume of simu- 
lation cell is 12 X 12 X 12 unit 3 (1 unit is 0.543 nm), which has 
13,824 atoms. 

There is finite size effect in calculating thermal conductivity when 
the simulation cell is not big enough 1819 . As shown in Fig. 2, we 
calculated thermal conductivity of isotopic nanoscale 3D phononic 
crystal with different size of simulation cell by equilibrium molecular 
dynamic (EMD) method at 1000 K. The period length is set as 
2 units and the mass ratio is set as 2. The calculated value of thermal 
conductivity converges when the side length of cubic simulation cell 
is larger than 10 units. To overcome the finite size effect on the 
calculated thermal conductivity, we use the side length of simulation 
cell as 12 units in the following simulations. 

For comparison with the thermal conductivity of pure 28 Si at 
1000 K, we also calculate its value as 50 ± 2 W/m-K (the dash dot 
line in Fig. 3(a)). Our result is comparable to Schelling et a/.'s results 
of MD simulation, 61 W/m-K 18 . However, MD results can not 
exactly coincide with the experimental value of 28 Si at 1000 K 20 , 
around 30 W/m-K, because of the inaccuracies of semi- empirical 
potentials and the impurity of the sample in measurements. This 
non-coincidence has little effect on the comparing MD results cal- 
culated using same potential parameters. 

We calculated the thermal conductivity of 3D phononic crystal 
with different period length, where the mass ratio R is 2. As shown in 
Fig. 3(a), the thermal conductivity rapidly decreases as the period 
length increases. The smallest value of thermal conductivity is 
2.14 W/m-K, which is only 4.3% of pure 28 Si calculated by EMD 
method. Simkin and Mahan show 12 that increasing period length 
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Figure 2 | Thermal conductivity versus the side length of simulation cell. 

The mass ratio of the 3D isotopic phononic crystal is 2 and the period 
length is 2 units. The error bars are calculated from 16 simulations with 
different initial conditions. All values are calculated at 1000 K which is 
larger than the Debye temperature, T D , of Si (—658 K). 

may increase the amount of band folding and decrease the average 
velocity in the superlattice, resulting in a decrease of thermal con- 
ductivity. The phonon mean free path of Si is around 60 nm at 
1000 K 18 , which is much longer than the period length in our simu- 
lation. The tendency of thermal conductivity in Fig. 3(a) is consistent 
with Simkin and Manhan's results when the phonon mean free path 
is larger than period length 12 . 

Another way to modulate the phonon transport in the 3D pho- 
nonic crystal is to vary the mass. The mass of impurity atoms could 
perturb the phonon density of state and phonon dispersion curves, 
which can affect the group velocities. 

Fig. 3(b) shows the dependence of thermal conductivity on the mass 
ratio, where the period length is 12 units. Our results indicate that 
thermal conductivity rapidly decreases as the mass ratio increases 
from 1 to 6. The heaviest Si isotope atoms produced is 43 Si 21 . where 
the mass ratio R corresponds to 1.5. The value of thermal conductivity 
is 4.2 W/m-K, that is, 8.4 percent of pure 28 Si (50 W/m-K). 



(a) 



(b) 





(c) 




(d) 




Figure 1 | The structures of the isotopic nanoscale 3D phononic crystals - 
three dimensional periodic arrangements of 28 Si and M Si atoms. From (a) 
to (d), the period lengths of those 3D phononic crystals are 2, 4, 6 and 
12 units, respectively. The lattice constant is 0.543 nm, that is, 1 unit 
represents 0.543 nm. In simulations, the periodic boundary condition is 
applied in all three directions. 
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Figure 3 | (a) Thermal conductivity versus the period length of isotopic 
nanoscale 3D phononic crystal of Si. The mass ratio is 2. The dash dot line 
corresponds to the molecular dynamic result of thermal conductivity of 
pure 28 Si. (b) Thermal conductivity versus the mass ratio of isotopic 
nanoscale 3D phononic crystals of Si. The period length is 12 units. All 
values are calculated at 1000 K. 
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Artificial Si isotopic atoms are used here to explore the mass 
influence on thermal transport and show the trend of large mass 
effects. The smallest value of thermal conductivity is 0.54 W/m-K 
(when R is 6), which is only 1.1% of pure 28 Si. Artificial M Si atoms can 
be looked as other atoms, such as 54 Fe 22 . When there are other kind 
atoms, the system is more complicated. That is, the mass is not the 
only factor involved. The bond strength and lattice relaxations must 
play a role, which is not studied in this letter. 

In a 28 Si 182 M Si 10 quasi-lD supercell with 10 M Si atoms (5%) ran- 
domly distributed, Gibbons and Estreicher 22 found the thermal con- 
ductivity decreased first and reached a minimum when the mass 
ratio was "two", and then the thermal conductivity increased as 
the increase of M. However, they stated that they cannot comment 
about the reasons for this minimum and do not know if the factor 
"two" remains valid for concentrations other than 5%. Different 
from randomly distributed, M Si in 3D phononic crystal is periodic 
distributed in 28 Si. The concentration of M Si (50%) is much bigger 
than 5%. Our results show that the thermal conductivity of 3D pho- 
nonic crystal decreases monotonously with increase of M. This is 
coincidence of the monotonous decrease of group velocities 
(Fig. 4(b)). 

As shown above, changing the mass of impurity atoms and the 
period lengths are two effective ways in modulating the thermal 
conductivity. To find the mechanism in the decrease of thermal 
conductivity of 3D phononic crystal, the phonon dispersion curves 
are calculated through classical lattice dynamics. We calculated the 
dispersion curves by general utility lattice program (GULP) 23 , and 
Stillinger- Weber potential 24 which is the same atom interaction as in 
our MD simulation. 




r a R 

Figure 4 | (a) Acoustic and partial optical branches along the [1, 1, 1] 
direction. Dispersion curves in red (black) correspond to the nanoscale 3D 
phononic crystal with 2 (4) unit cells in period length. The mass ratios are 
2. It is clearly shown close to R point that the group velocities decrease as 
the period length increases, which causes the reduction of the thermal 
conductivity, (b) Acoustic branches along the [1, 1, 1] direction. The mass 
ratio of 3D phononic crystals changes from 1 to 2.5. Different color is 
referred to different mass ratio. The period lengths of 3D phononic crystals 
are kept same, 2 units. The dispersion curves are affected by the mass of 
impurity atoms and the group velocities decrease as the mass of impurity 
atoms increase, which contributes to the decrease of thermal conductivity. 



Fig. 4(a) shows acoustic branches and partial optical branches of 
the dispersion curves of the 3D phononic crystal, where period 
lengths are different, and the mass ratio is the same. The Brillouin 
zone with 2 units in period length is eight times as large as that with 
4 units in period length. In calculating the dispersion curves of struc- 
ture with 2 units in period length, we use a larger unit cell whose size 
is as the same as that of structure with 4 units in period length. 

As the optical phonons contribute less to the thermal conductivity 
due to the lower group velocities, we focused on the acoustic pho- 
nons. It is clearly shown in Fig. 4(a) close to R point that the group 
velocities decrease as the period length increase, which causes the 
reduction of the thermal conductivity (shown in Fig. 3(a)). Fig. 4(b) 
shows acoustic branches of the 3D phononic crystal with different 
mass ratios, where the period length is 2 units. The dispersion curves 
are affected by the mass of impurity atoms and the group velocities 
decrease as the mass of impurity atoms increase, which contributes to 
the decrease of thermal conductivity (shown in Fig. 3(b)). 

To understand more about the underlying physical mechanism of 
thermal conductivity reduction, we carry out a vibration eigen-mode 
analysis on 3D phononic crystals. The mode localization can be 
qualitatively characterized by the normalized inverse participation 
ratio (NIPR) 25 . The NIPR for phonon mode k is defined through the 
normalized eigenvector u k 

2=1 \ a =l / 
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Figure 5 | The normalized inverse participation ratio (NIPR) spectra. 

NIPR is calculated based on Eq. (1). The larger of the value of NIPR the 
more localized of a phonon mode, (a) NIPR spectra of 3D phononic 
crystals with different period length. The left and right panels are 
corresponding to 3D phononic crystals with 2 unit and 4 unit period 
length, respectively. The mass ratios are the same as 2. (b) NIPR spectra of 
3D phononic crystals with different mass ratio, R. The period lengths are 
the same as 2 units. The upper left panel (R = 1) corresponds to pure Si. 
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where N is the total number of atoms. When there are less atoms 
participating in the motion, the phonon mode has a larger NIPR 
value. For example, NIPR is N when there is only one atom vibrates 
in the localized mode. When all atoms participate in the motion, 
NIPR is calculated out as 1. That is, the larger of the value of NIPR 
the more localized of a phonon mode. 

Fig. 5(a) shows the NIPR spectra of 3D phononic crystal with 
two different period lengths, where the mass ratio R is 2. 
Obviously, the NIPR for 3D phononic crystal with 4 units in per- 
iod length has larger values than that with 2 units in period length. 
That is, there are more modes localized in 3D phononic crystal 
with 4 units in period length, which also leads to a reduction of its 
thermal conductivity. Fig. 5(b) shows the NIPR spectra of 3D 
phononic crystal with different mass ratio R, where the period 
length is 2 units. The values of NIPR for R = 1.5, R = 2 and R 
= 2.5 are close to each other and larger than pure Si (R = 1), 
which show that the isotopic atoms could cause more localizations. 
In Fig. 5(b), band gaps appear in the high frequency part 
(>12 THz) of spectra for R = 2 and R = 2.5, and the band gaps 
for R = 2.5 are wider than those for R = 2. That is, the phonons 
with frequency in the range of band gaps cannot exist in the 3D 
phononic crystal. When R is smaller than 1.5 (corresponding to 
43 Si), band gaps are not observed in 3D isotopic phononic crystal. 
So, band gaps are minor effects when the mass difference is not big 
enough. 

Besides decreasing group velocity and phonon modes localization, 
interface disorder is another factor in reducing the thermal conduc- 
tivity for superlattices 11 and phononic crystal 16 . In our simulation, 
28 Si and M Si atoms has a perfect interface which do not include the 
disorder effect. There will be a further deduction in thermal conduc- 
tivity after considering the interface disorder. So, the disorder effect 
will be studied in the future work. 

Discussion 

Using the Green-Kubo method, we have calculated thermal conduct- 
ivities of isotopic nanoscale 3D phononic crystals, where 28 Si and M Si 
atoms are assembled periodically in the three directions. Results 
show that the thermal conductivity decreases as the increasing of 
period length from 1 nm to 6 nm. The thermal conductivity of struc- 
ture with 6 nm period length and 2 in mass ratio is 2.14 W/m-K at 
1000 K, which is only 4.3% of pure 28 Si and can lead to a larger ZT 
than unity. 

Moreover, the thermal conductivity rapidly decreases as the 
mass ratio increases. The phonon localizations and bandgaps at 
high frequency in the 3D phononic crystal are shown clearly in 
the spectra of the normalized inverse participation ratio. The 
appearance of band gaps blocks a range of frequency of phonon 
modes and flattens phonon dispersion curves. The phonon dis- 
persion curves show the phonon group velocities decrease in 3D 
phononic crystal. In a word, the decrease of thermal conductivity 
in 3D phononic crystal is attributed to both the decrease of group 
velocities and the localization. 

Thermal conductivity is mainly contributed from acoustic pho- 
nons. If the bandgaps of the 3D phononic crystals can exist at low 
frequencies, there will be greater reduction of the thermal conduc- 
tivity. However, manipulating band gaps to low frequencies in 
nanoscale 3D phononic crystal is a challenging work which is worthy 
of effort in the future. There are advances in obtaining the nanoscale 
2D phononic crystal. However, the limitation in fabricating nanos- 
cale 3D phononic crystal is still challenging nowadays. 

Methods 

The Green-Kubo method, equilibrium molecular dynamics (MD), is employed to 
calculate the thermal conductivities of 3D phononic crystal at 1000 K. MD simulation 
is a popular method in calculating thermal conductivity at high temperature 18,26 ' 27 . In 
this letter, we focus on the thermal conductivity of 3D phononic crystal at 1000 K, 
which is larger than the Debye temperature, T D , of Si (—658 K) 28 . 



In simulations, the periodic boundary condition is applied in all three directions. 
To derive the force term, we use Stillinger-Weber (SW) potential for Si 24 , which 
includes both two-body and three-body potential terms. The SW potential has been 
used widely to study the thermal properties of Si bulk material 14 ' 26 ' 29 for its accurate fit 
for experimental results on the thermal expansion coefficients. The heat current is 
defined as 18 

7,(o = E *<« + 1 E 7 -) M + E h iHim) (2) 

i l} i¥=j ijk 

where Fy and F ^denote the two-body and three-body force, respectively. Thermal 
conductivity is calculated from the Green-Kubo formula 30 

K= 3k^v[ <m ^ (3) 

where k is thermal conductivity, ks is the Boltzmann constant, V is the system 
volume, T is the temperature, and the angular bracket denotes an ensemble average. 

Generally, the temperature in MD simulation, T MD , is calculated from the kinetic 
energy of atoms according to the Boltzmann distribution: 

n 1 3 

{E)=Y J ~mv 2 i = -Nk B T MD (4) 
i=\ l 1 

where (E) is the total kinetic energy, v t is the velocity, m is the atomic mass, N is the 
number of particles in the system, and k B is the Boltzmann constant. This equation is 
valid at high temperature (T^>T D , T D is the Debye temperature). 

Numerically, velocity Verlet algorithm is employed to integrate equations of 
motion, and each MD step is set as 1.0 fs. Firstly, canonical ensemble MD with 
langevin heat reservoir runs for 2 20 steps to equilibrate the whole system at 1000 K. 
Then, microcanonical ensemble (NVE) MD runs for another 2 25 steps (33.5 ns). 
Meanwhile, heat current is recorded at each step. At the end, the thermal conductivity 
is calculated by Eq. (3). In the calculation of thermal conductivity, the integration is 
from zero to a cut-off time which is determined by "first avalanche" method 31 . The 
final result is averaged over sixteen realizations with different initial conditions to 
satisfy ergodicity. 
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